Draft version June 16, 2009 

Preprint typeset using I£T^X style cmulateapj v. 10/09/06 



O 

o 

(N 



AN ESTIMATE OF THE PRIMORDIAL NON-GAUSSIANITY PARAMETER i NL USING THE NEEDLET 

BISPECTRUM FROM WMAP 

0YSTEIN RUDJORD 1 ' 2 , FRODE K. HANSEN 2 , XlAOHONG LAN 3 , 
MlCHELE LlGUORI 4 , DOMENICO MARINUCCI 3 , SABINO MATARRESE 5 
Draft version June 16, 2009 

ABSTRACT 

We use the full bispectrum of spherical needlets applied to the WMAP data of the cosmic microwave 
background as an estimator for the primordial non-Gaussianity parameter Jnl- We use needlet scales 
up to ^ max = 1000 and the KQ75 galactic cut and find Jnl = 84 ± 40 corrected for point source 
bias. We also introduce a set of consistency tests to validate our results against the possible influence 
of foreground residuals or systematic errors. In particular, fluctuations in the value of Jnl obtained 
from different frequency channels, different masks and different multipoles are tested against simulated 
maps. All variations in /jvl estimates are found statistically consistent with simulations. 
Subject headings: cosmic microwave background — cosmology: observations — methods: statistical 
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1. INTRODUCTION 

The theory of inflation predicts the fluctuations in the 
Cosmic Microwave Background (CMB) to be close to 
Gaussian distributed. However, a small degree of non- 
Gaussianity is generally present in all the inflationary 
scenarios. The primordial non- Gaussian signal predicted 
by many models can be parametrized in the form: 

$(x) = * L (x) + /nl - (*i(x)» , (1) 

where $(x) is the primordial curvature perturbation 
field at the end of inflation and $d(x) is the Gaus- 
sian part of the perturbation. The dimensionless pa- 
rameter /nl describes the amplitude of non-Gaussianity. 
The non-Gaussian part of the primordial curvature 
perturbation is a local functional of the Gaussian 
part and for this reason this kind of parametriza- 
tion is often referred to as local non-Gaussianity. Lo- 
cal non-Gaussianity is predicted to arise from stan- 
dard single-field s low-ro ll inflation (|Acauaviva et al.~l 
120031 : iMaldacena I l2003h as well as from alternative 
inflationary scenarios for the gener ation of primor- 
dial perturbations, like the curvaton (lEnqvist fc Sloth] 
12002 iLvth fc Wands 1 12001 iMoroi fc Takahashi I l200lT 
or in homogeneous (pre] reheating models (|Dvali et al. I 
120041: iKolb et al. 1 120051 120061 ). or even from alterna- 
tives to inflation, such as ekpyrotic and cyclic models 
(jKovama et al. 1 120071 : iBuchbinder et al. 1 12008^ Other 
models, such as DBI inflation ( Alishahi ha et al. Il2004h 
and ghost inflation (|Arkani-Hamed et al. 1 120041 ), pre- 
dict a different kind of primordial non-Gaussianity, 
called " equilateral" , because the three point function 
for this kind of non-Gaussianity is peaked on equilat- 
eral configurations, in which the lengths of the three 
wavevectors forming a triangle in Fourier space are equal 

1 email: oystein.rudjord@astro.uio.no 

2 Institute of Theoretical Astrophysics, University of Oslo, P.O. 
Box 1029 Blindern, N-0315 Oslo, Norway 

3 Dipartimento di Matematica, Universita di Roma 'Tor Ver- 
gata', Via della Ricerca Scientifica 1, 1-00133 Roma, Italy 

4 Department of Applied Mathematics and Theoretical Physics, 
Centre for Mathematical Sciences, University of Cambridge, 
Wilberfoce Road, Cambridge, CB3 0WA, United Kingdom 

5 Dipartimento di Fisica, G. Galilei, Universita di Padova and 
INFN, Sezione di Padova, via Marzolo 8,1-35131 Padova, Italy 



(|Creminelli et al.~l 120061 ) . In this paper we will focus 
only on non-Gaussianity of the local type, described by 
equation [TJ The interesting aspect of primordial non- 
Gaussianity is that the expected non-Gaussian amplitude 
/nl varies significantly from model to model. Putting 
experimental bounds on /nl is then equivalent to con- 
straining primordial scenarios of inflation. For example 
standard single-field slow-roll i nflation predicts /nl ~ 
10~ 2 at the end of inflation (|Acauaviva et al. 1 120031 : 
lMaldacena1l2003l ) (and therefore a final value ~ unity af- 
ter general relativistic s econd-order perturbat ion effects 
are taken into account (Ba rtolo et al. ll2b04bl fcT)). Such 
a small value is not experimentally detectable and for 
this reason an eventual detection of a Gaussian signal in 
present and forthcoming CMB data will rule out single- 
field slow-roll inflation as a viable scenario. Motivated 
by these considerations many groups have attempted to 
measure /nl using CMB datasets, and WMAP data in 
particular. 

A detection of n on-zero /jyx at more th an the 2a 
level was found by (|Yadav fc Wandelt 1 120081 using the 
WMAP data with the KpO galactic cut. The WMAP 
team found similar values but stating that only the 
value obtained with the slightly larger KQ75 galactic 
cut is reliable due to possible foreground residuals. In 
this case a value oi /nl = 51 ± 3 2 was found. In 
both these cases, an extended version (|Creminelli et al. 1 
l27)fflljYadav fc Wandelt I l2008t lYadav et al. I I2007H) of 
the (jKomatsu et al. Il2005h (KSW-method) based on the 
full bispectrum was us ed. C onsistent results were fo und 
by (jCurto et al. Il2008f) and (|Pietrobon et al. Il2008a[) us- 
ing parts of the bispectrum of spherica l mexican hat 
wavelets ([Martinez- Gonzalez et al.ni2002[ ) and the skew- 
ness of needle t coefficients. A re cent estimate has now 
been made by ([Smith et al. II2008T ) obtaining the smallest 
error bars on Jnl so far findi ng f^L = 38 ± 21 

Recently it was shown in (|Lan fc MarinuccTl l2008al) 
that needlet coefficients can be used to construct statis- 
tics which are directly related to the bispectrum. These 
statistics share most of the useful properties of the bis- 
pectrum, while at the same time they do present im- 
portant advantages, especially in terms of robustness 
to masked data and computational rapidity. Motivated 
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from these results, in this paper we will use the full bis- 
pectrum of needlet coefficients as an estimator of Jnl, 
introducing moreover a set of consistency tests to check 
the stability of our findings. In particular, we shall in- 
vestigate whether changes in the estimated values of Jn l 
using different galactic cuts, different frequency bands 
and different multipoles are within the variations ex- 
pected from simulations. Of course, variations in Jnl 
among these different cases which are significantly larger 
than statistical fluctuations might point out the pres- 
ence of foreground residuals or other systematic effects 
that could have biased the estimate of Jnl- 

The plan of the paper is as follows. In section [21 we 
describe the data used in the analysis. Then in section 
[3l the needlets and the needlet bispectrum estimator are 
described in detail. Finally, the results on the WMAP 
data are presented in section [5] and conclusions are made 
in section 

2. DATA 

For this analysis we used the noise weighted average of 
the V and W frequency band of the WMAP 5 year CMB 
map, as well as the corresponding instrumental beam and 
noise properties. We have also performed the analysis on 
the individual Q (41 GHz), V (61 GHz) and W (94 GHz) 
bands. For masking out galactic foregrounds and point 
sources, we used the KQ75 and KQ85 mask supplied by 
the WMAP team. For particular cases, we also used the 
much smaller Kpl2 mask (maintaining 94% of the sky) 
as well as an extended KQ75+ mask. The KQ75+ mask 
is constructed from the KQ75 mask, extending the mask 
with 5 degrees along the rim, maintaining a total of 63% 
of the sky. We have used the maps at Healpix 6 resolution 
N s ide = 512. 

3. METHOD 

3.1. Spherical needlets 

Needlets are a new form of (second generation) spheri- 
cal wavel ets, which were introduc ed into functional anal- 
ysis by ( (|Narcowich et al.ll2006al lbT)) and have attracted 
a lot of attention in the cosmological literature here- 
after. The possibility to use needlets for the statis- 
tical analysis of spherical random fields, with a view 
to C MB applications, is first discussed in (jBaldi et al.l 
2006), where the stochastic properties of needlet co- 
efficients are established and their possible roles for 
data analysis (spectrum estimation, Gaussianity testing) 
are described; further mathemat ical properties where 
then given in (|Baldi et al. I I2007D . The first applica- 
tion to CMB data, in particular, for the analysis of 
cross-correlation of CMB and Large Scale Structure data 
was provided by (|Pietrobon et al. 11200 61: a general pre- 
sen tation of the method fo r a CMB a udience is given 
in rtMarinucci et al. I l2008f ), while in (jGuilloux et al. I 
l2007fl the properties of different weighting schemes are 
investigated and compared. Further papers have ap- 
plied needlets on CMB data, for issues such as map- 
making, spectru m estimation, detection of features 
and anisotropics (([Fay et al. 1120 08b: Dcl abrouille et ai~l 
120081: iFav et al. ll2008aHPietrobon et al. Il2008bf) 1: more 
recently, needlets have also been considered for the 

6 http://healpix.jpl.nasa.gov 



analysis of directi onal data, with a view to high en- 
ergy cosmic rays ( (|Baldi et al. 1120081)1 and for the anal- 
ysis of polarizatio n data ( (|Geller fc Marinucci I 120081 : 
iGeller et al. 1 1200811. whereas extensio ns to the so-called 
Mexic an needlets case are discussed by ()Geller fe Mavelil 
l2007blR), their stochastic properties bei ng established 



in (|Lan fe Marinucci 



nastic propertie s pern 
1 l20081ilM^eTi1[200l . 



The spherical needlet (function) is defined as 



E 6 (^) E Y UWU7k) ; (2) 



here, 7 is a direction on the sphere, and j is the frequency 
(multipole range) of the needlet and Xjk is a normalizing 
factor. The points {7^} can be identified with the pixel 
centres in the HealPix pixelization scheme. The number 
B defines the needlet basis such that only multipoles in 
the range i = [B^~ l , _B J+1 ] are included, i.e. the function 
b(£/Bi) is zero outside this range. For details in the func- 
tiona l form of b^/B^). please refer to (|Marinucci et al.l 
2008!) and references therein. 

The advantages of needlets have already been discussed 
in several papers in the literature; in short, we recall that 
needlets do not rely on any tangent plane approximation; 
they are computationally very convenient, and inherently 
adapted to standard packages such as HealPix; they al- 
low for a simple reconstruction formula (a property which 
is not shared by other wavelets systems) ; they are quasi- 
exponentially (i.e. faster than any polynomial) concen- 
trated in pixel space. Moreover, needlets are exactly lo- 
calized on a finite number of multipoles (the width of 
this support is explicitly known and can be specified as 
an input parameter, see Eq. [2])). 

Random needlet coefficients can be shown to be asymp- 
totically uncorrelated (and hence, in the Gaussian case, 
independent) at any fixed angular distance, when the fre- 
quency increases. This capital proper ty is in general not 
shared by other wave let systems (see ( (|Lan fc Marinuccil 
2008b; M avelil [20081 )11 and can be exploited in several 
statistical procedures, as it allows one to treat needlet 
coefficients as a sample of independent and identically 
distributed coefficients on small scales, at least under the 
Gaussianity assumption. 

In this paper, for notational simplicity we shall take 
random needlet coefficients to be 

I 1 

Pjk = E & (#j) E a e™Ye m (%) = ^b e a em Y em (%) 

£ m—~£ ira 

Here j denotes the frequency of the coefficient and 7^ 
denotes the direction on the sky. We are dropping a nor- 
malizing term yj\jk\ this comes at no cost, because in 
this paper needlet coefficients always appear after nor- 
malization with their own standard deviation, so that 
this deterministic factor cancels. From the notational 
point of view, however, this allows a major simplifica- 
tion, as it permits to avoid considering different weights 
at different frequencies j. As before, the index k can in 
practice be the pixel number on the HEALPix grid. 

3.2. The needlets bispectrum 

Starting from some highly influential papers 
( (|Hu et al. I [2OOII iKomatsu fc Spergell [200ll) 1. the 
bispectrum has emerged in the last decade as the most 
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promising statistics for the detection of non-Gaussianity 
in CMB data. To fix notation, we recall that, under the 
assumption of statistical isot ropy for CMB radiation, we 
must have ( l|Hu et al. 1120011 )) 



\Q|i?7li ^timi ^£37713 / 



ix li h 

m\ m<i 



D 



where on the right hand side we have introduced 
the Wigner's 3j coefficients, which are different from 
zero only for configurations of h^hi h which s atisfy 
the triangle conditions (see again (IH11 et al. I 120011 : 
iKomatsu fe Spergel I l200lh ) . The unreduced bispectrum 
£?^ 2 £ 3 is identically zero in the Gaussian case; under 
non-Gaussianity, it can be estimated by 



B, 



E 



7711777.27713 



k h 4 

mi ?Ti2 7713 



It was shown by ( (IKomatsu fe Spergel I l200l[ )) that, in 
the idealistic circumstance with the absence of masked 
regions, the bispectrum can constrain non-Gaussianity 
very efficiently, with a signal-to-noise ratio equal to 
unity for Jnl smaller than 10 at the Planck resolu- 
tion. In the presence of masked regions, however, these 
properties deteriorate consistently; many statistical so- 
lutions have been d iscussed so far, see for instance 
(lYadav et al. I 120071 : lYadav fc Wandelt I HH) for the 
most recent developments. A large literature has also 
focussed on the determinations of the multipole configu- 
rations where the signal-to-noise ratio should be expected 
to be stronger, in view of a given model: see for instance 
dBabich et al. ll2004ICabella et al. ll2006HBartolo et al. I 
I2004al iMarinucci I l2006t lYadav fe Wandelt 1 120081) and 
many others. 

Our purpose in this paper is to combine ideas from 
the bispectrum and the needlets literature, to propose a 
needlet bispectrum method to test non-Gaussianity and 
estimate the nonlinearity parameter /jvl- More precisely, 
we suggest to focus on the needlet bispectrum, defined 

by 



eJ jik <J j2kVj 3 k 



where crj = w < fi jk > is the standard deviation of 

Pjk- The needlet bispectr um was first considered in 
( ()Lan fe Marinuccil l2008 a)). and we refer to that pa- 
per for more discussion and details on its mathemati- 
cal properties; the use o f needlets for a non-Gau ssianity 
test is also proposed in (iPietrobon et al. Il2008al ). where 
the focus is instead on the skewness of the coefficients 
(which can be viewed as a special case of the bispectrum, 
for ji = ]2 = J3). Of course, many other papers had 
previously used wavelet-related techniques to search for 
non-Gaussianity in CMB, see for instance (IVielva et al. I 
[200llCruz et al. II2007L l200l iMcEwen et al. Il2006l) ~ 

In short, to understand the motivations of our pro- 
posals, note that, denoting by N the cardinality of the 
points k (i.e., the number of points in the pixclization 
scheme, so that 4n/N provides an approximation for the 
pixel area), and neglecting for simplicity beam factors, 
we have approximately 

-^^203lkf3]2kf3hk = 



iimi £2^1 £3^3 

^17711 ^£2^2 %"i3 

Y ixmi (%)Yt 2m2 (%)Yi 3m3 (%) 



~E K^)H^)H-j^;)ai imi ae2 m 2ae 3m3 



Write 



, tx h h\ /(2^i + l)(2£ 2 + l)(24 + l) 

^1^2^3 = I " 



4tt 



hence wc obtain 



E ^ OJ! )K DJ 2 )K Dj 3 ) 

iirrii 

( h 4 is 1 , 

x ae imi ae2ni2 a e 3 rn 3 I mi m2 m ^ ) ne^ts 



Eh 



From the previous computations, it should be clear that 
the needlets bispectrum can be viewed as a smoothed 
and normalized form of the standard bispectrum estima- 
tor. As usual with wavelet techniques, the advantage is 
that, while the standard bispectrum is known to be heav- 
ily affected by the presence of masked regions, needlet 
coefficients are much more robust and consequently the 
needlet bispectrum makes up a more reliable statistics 
even for incomplete maps. Furthermore, the needlet bis- 
pectrum is computationally very convenient, as it does 
not require the evaluation of Wigner's 3j coefficients, 
which is extremely time-consuming. 

From the mathematical point of view, further prop- 
erties of the needlets b ispectrum are discussed by 
(|Lan fc Marinuccil l2008af) : in particular, it is shown 

that, after normalization, Ij 1 j 2 j 3 is asymptotically Gaus- 
sian as the frequency increases. Furthermore, it can 
be shown that (under idealistic experimental circum- 
stances) the values of the needlet bispectrum are asymp- 
totically independent over different frequencies, so that 
chi-square statistics can be suitably implemented. In 
()Lan fc Marinuccil l2008a( ) . some analytic discussion on 
the power properties of the needlet bispectrum for a pure 
Sachs- Wolfe model were also provided, showing that its 
expected valued diverges to infinity at high frequencies. 
Although those results were derived in a mathematical 
setting and did not take into account many features of 
CMB data, the simulations in the present paper show (in 
a much more realistic setting) that this procedure does 
have very satisfactory power properties in the presence 
of non-Gaussianity. 
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3-3. Jnl estimator 

We will now use the needlets bispectrum for estimating 
Jnl by a x 2 minimization procedure. We define X 2 (Jnl) 
as 

X 2 (fN L ) = d T (f NL )C- 1 d(f NL ), 

where the elements d{ of the data vector d are defined as 
di = Ij 1 j 2 j 3 (observed) — (Ij 1 j 2 j 3 )(fNL) for all combina- 
tions of [ji, j2, ja) satisfying the triangle condition. Here 
Ij 1 j 2 j 3 (observed) is the needlets bispectrum of the ob- 
served data and {Ij 1 j 2 j 3 )(fNL) is the expectation value 
of the needlet bispectrum for a given value of /at l ■ The 
correlation matrix C is given by 

dj = {didj) - (dijidj). 

The correlation matrix is obtained from Gaussian sim- 
ulations. In order to avoid cumbersome numeric grid- 
calculations of {Ij 1 j 2 j 3 } {/nl), we seek an expression with 
a more explicit dependency of /nl- In order to arrive at 
such an expression, we write out again the needlets bis- 
pectrum as 



npix 



/3jikf3j 2 kf3j 3 ii 



VjlkO-j 2 kO- j3 k 



V 



k e 1 £ 2 £ 3mi m 2m3 

x Y iimi (%)Y i2m2 (%)Yt sm3 (%) 



ei mi t - t, t2m2 ^t^m^ 



(3) 

As usual, the non-Gaussian aim's are assumed to be a 
combination of a linear (Gaussian) and a non-linear term: 
o-im = af m + fNLa-i'm- This allows us to write the three- 
point correlations in a^ m 's as 

(^limi &i 2 rn 2 &i 3 m 3 ) 

— / n G n G n G \-Uf ( / n NG n G n G \ 

— \ a iim 1 a £ 2 m 2 a £ 3 m 3 / + J NL { \ a l imi a l 2 m 2 a i 3 m 3 / 
+ ( a f 1 m 1 a i' 2 m 2 a f 3 m 3 ) + { a ? imi a f 2 m 2 a f 3 m 3 ) ) 



+ o((«L G ) 2 ) 



(4) 



The non-linear terms are assumed to be small, and thus 
we will neglect the higher order terms, 0{(af^) 2 ) rj 0. 
We will also neglect the pure Gaussian term, since the 
three point correlation function of a Gaussian field is 
zero. We insert the remaining terms into eq. [3] and find 
the mean value: 



( I 3U2j 3 ) (Jnl) - Jnl | 



I npix R G 0NG0G 
PjikPj 2 k Pjsk 

a jlkO-j 2 kO-j 3 k 



/npix R NG oG oG 
\ "* Pjik Pj 2 kPj 3 k 

\u a hk°~32kO-j 3 k 



/npix R G oG oNG 
P hkPj 2 kPj 3 k 



o-j 1 kO- 32 kO- 33 k 



= Jnl { Ij x j 2 j 3 



NG,G,G 
3233 



jG,NG,G 
hh33 



where we have defined the quantity 



g,g,ng\ 

113233 I 



(5) 



rNG,G,G 
3ih33 



+ I 



jG,NG,G 
"313233 



G,G,NG 
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(dashed) plotted along j\ = j while j% = 25 and 
J3 = 25, compared to the average bispectrum from 300 simulations 
with /jvl = 500 (full line) and average bispectrum from 10000 
Gaussian simulations (dotted). 

which does not depend on f^L to the first order in a^ . 
Figure Q] shows a plot of a average needlet bispectrum 
Cii ja h ) (Inl) from 300 simulations with /jvl = 500. 
Here the bispectrum is plotted along one of the indices 
3i = i> while the two other indices ]2 — 25 and j'3 = 25 
are kept constant. In the same plot is also our first order 
approximation, 500 x (Ij 1 j 2 j 3 )- As we see, this approxi- 
mation is fairly g ood for Jnl = 500, a nd based on pre- 
vious estimates ( (jKomatsu et al.ll2008l )). we will assume 
that Jnl does not have a value significantly higher than 
this. 

We can now write the elements of the data vector d of 



the x as di = Ij 



- Jnl {I jt 3233 ) ■ In order to estimate 



Jnl , we need to find the value of Jnl that minimizes the 



giving 



d X 2 (fNL) 



df. 



0. 



NL 



(G) 



[hihh) c ll hhh (observed) 
fNL = ; — ; ; • (7) 



Ijihh ) ^ 1 ( I313233 



3.4. 



The procedure to estimate Jnl 

We will now show the full procedure we have used to 
estimate Jnl using the needlet bispectrum. 

1. Generate 10000 simulations of Gaussian sky maps 
using the best fit WMAP 5 year power spectrum. 
These are smoothed with an instrumental beam 
and noise is added. The maps arc also multiplied 
with a mask for galactic cut, in order to remove 
foregrounds. A needlet transform is applied and 
the standard deviation ajk of the needlet coeffi- 
cients f3jk are found. This standard deviation is 
needed to find the needlet bispectrum as seen from 
eq. ©. 

2. Produce another 120000 Gaussian simulations. 
Mask, beam and noise properties are applied as 
above. The needlet coefficients are found and used 
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Fig. 2. — A histogram (full line) of /jvz, estimates of 10000 
Gaussian simulations of the V+W channel plotted together with 
a Gaussian fit (dashed). As expected the /jvl values form an 
approximate Gaussian distribution with a mean value of /jvl = 0. 

to obtain the needlet bispectra, Ij 1 j 2 j 3 ■ These bis- 
pectra are used to find the covariance matrix C. 
This converges very slowly, thus the need for such 
a large number of simulations. 

3. Generate 300 non-Gaussian simulations 

((Lig uori et al. I [2007)) to find the mean first 

order non-Gaussian bispectrum, (Ij 1 j 2 j s )- 

4. Obtain the needlet bispectrum from the data and 
estimate /jvz, using eq. [7] 

5. Generate a set of 10000 simulated Gaussian maps 
and estimate Jnl in the same manner from these 
maps in order to obtain the error bars on Jnl- 
This set of estimated /jvl values form a Gaussian 
distribution around /jvl = 0. Figure [5] shows a 
histogram of the 10000 Jnl of the V+W frequency 
channel estimates plotted together with a Gaussian 
distribution. We see that the distribution of /nl 
estimates is very close to Gaussian, so we will use 
the standard deviation as a measure of the uncer- 
tainty of the estimate. 

4. RESULTS 
4.1. Estimates of f^L 

We used the above procedure to estimate Jnl from 
the WMAP data. The co-added V+W map as well as 
the single frequency bands Q, V and W were used. The 
estimated values of /jvl are listed in table Q] together 
with the lcr error bars found from the simulations. For 
the analysis we used both the KQ85 mask and the more 
aggressive KQ75 mask for galactic cut in order to study 
potential effects from residual foregrounds. 

We have used multipoles up to £ max = 1024 in the 
analysis. We tested different values of B in order to find 
the number of needlet scales from £ = 2 to £ = 1000 which 
would yield an invertible covariance matrix. We found 



that the maximum number of scales which could be used 
was 31 scales, using B = 1.2050. For comparison with 
the WMAP results we also used needlet scales including 
multipoles up to £ max = 500 and £ max = 700. In table 
[T] as well as in the following text, we will use ^ max = 
1000, 700, 500 to specify the highest multipole included 
in the highest needlet scale. Note that this number may 
differ slightly from 1000, 700 and 500 depending on the 
exact value of B specified. For £ max = 500, we used 30 
scales with B = 1.1828 and for ^ max = 700 we used 31 
scales with B = 1.1880. 

We see from the table that the best results obtained on 
the combined V+W band yields Jnl = 89 + 39 for the 
KQ75 cut and Jnl — 117 ± 36 using the smaller KQ85 
cut. We run simulations of unresolve d point sources 
based on the procedure described in (|Komatsu et al.l 
|2008|) and obtained a point source bias of A/at l = 5 ± 1 
for KQ75 and A/jy l = 7 ± 1 for KQ85 giving corrected 
values of f NL = 84 ± 40 and f NL = 110 ± 37. We see 
that even for KQ75 a zero value for /jvl is excluded at 
about 2a. In order to make sure that foregrounds are 
not influencing our results significantly, we also make an 
estimate on the much larger KQ75+ mask and obtain 
Jnl = 103 ± 41 or Jnl = 97 + 42 taking into account 
unresolved point sources. 



freq. channel 


mask 


@-max 


n i 


Jnl 


V+W 


KQ85 


700 


31 


156 ± 45 


V+W 


KQ75 


700 


32 


88 ±48 


V+W 


Kpl2 


1000 


31 


160 ± 30 


V+W 


KQ85 


1000 


31 


117 + 36 


V+W 


KQ75 


1000 


31 


89 ± 39 


V+W 


KQ75+ 


1000 


31 


103 ± 41 


V+W (Raw) 


KQ85 


1000 


31 


105 ± 36 


V+W (Raw) 


KQ75 


1000 


31 


83 ± 39 


V+W (Raw) 


KQ75+ 


1000 


31 


87 + 41 


V 


KQ75 


500 


30 


78 + 57 


V 


KQ85 


1000 


31 


100 ± 39 


V 


KQ75 


1000 


31 


105 ± 42 


V (Raw) 


KQ85 


1000 


31 


88 ± 39 


V (Raw) 


KQ75 


1000 


31 


100 ± 42 


W 


KQ75 


500 


30 


57 + 59 


w 


KQ85 


1000 


31 


79 + 42 


w 


KQ75 


1000 


31 


54 + 45 


W (Raw) 


KQ85 


1000 


31 


57 + 42 


W (Raw) 


KQ75 


1000 


31 


41 + 45 


Q 


KQ75 


500 


30 


47 + 59 


Q 


KQ85 


1000 


31 


33 + 42 


Q 


KQ75 


1000 


31 


9 + 44 


Q (Raw) 


KQ85 


1000 


31 


-64 ± 42 


Q (Raw) 


KQ75 


1000 


31 


-21 ± 44 


Combined V and W 


KQ75 


1000 


30 


76 ± 38 



TABLE 1 

The estimated values for /jvz, together with the la error 

BARS. 

As we see there is a large improvement in error bars 
from £ ma x = 700 to £ max = 1000. This may seem at first 
sight surprising taking into account the fact that this 
range is noise dominated. However, this result is not un- 
expected if we take into account the way the needlets 
are constructed. Indeed, £ maK = 700 means that no 
needlet scale using information above ^ max = 700 is in- 
cluded. Nevertheless, from the previous description of 
the needlet systems it is easy to see that the information 
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from multipoles close to the boundary value ^ max = 700 
will receive very little weight in general. Of course, the 
next needlet scale will contain information below as well 
as above ^ max = 700. Therefore, when extending the 
analysis to higher £'s we do not only exploit 300 more 
multipoles, but we are also able to extract better infor- 
mation from the multipoles below I = 700. 

Another test was performed to take advantage of the 
fact that the CMB should be independent of frequency, 
while the noise differs between the channels. A data vec- 
tor was composed from the needlet bispectrum of both 
the individual V and W frequency channels. 



IY. . 



tW 

jihte 



(8) 



The full covariance matrix in this case contains infor- 
mation about correlations between the frequencies, and 
should therefore enable us to get smaller error-bars on 
/jvl- However, for this analysis it was necessary to use 
only 30 needlet scales from each frequency channel in 
order to get an invertible covariance matrix. And the 
result (shown in the bottom row of table [TJ) was not a 
large improvement from the analysis of the VW band at 
31 needlet scales. However this is our estimate for /jvl 
with the smallest error-bars while using the KQ75 mask. 

For the B = 1.2050 case for the V+W band with the 
KQ75 mask, we have also investigated the change in /jvl 
as a function of the number of needlet scales included. 
We thus included only the first 25 needlet scales (up to 
< max = 324), then the 26 first scales (up to £ max = 390 
and so on up to all 31 scales. The results are presented in 
table [21 As expected, we see that the error bars are de- 
creasing with increasing £ max . Differently from our case, 
in the optimal bispectrum estimation performed by the 
WMAP team and other groups error bars saturates ear- 
lier than lmax = 1000 because the full inverse covariance 
weighting scheme is not implemented and an approxi- 
mation is used (whereas in this case the Monte Carlo 
approach used to estimate the bispectrum automatically 
accounts for this issue). Note however the WMAP error 
bars at lmax = 700 are still smaller than ours at lmax = 
1000 because we don't implement a minimum variance 
estimator and thus we don't saturate the Rao-Cramer 
bound. Moreover an optimal bispectrum estimator with 
full inverse covariance weighting has been very recently 
implemented by Smith et al. 2008. 

4.2. Consistency checks 

We see from these results that the estimates using the 
KQ85 mask differs notably from the estimates using the 
KQ75 mask. This is particularly the case for the V + W 
channel, when only considering scales up to £ m ax = 700. 
This estimate when using the KQ85 mask (/jvl = 156) 
is much higher than the estimate found from the same 
map, using the KQ75 mask (/jvl = 88). We are there- 
fore motivated to study simulations to find how often 
a change of mask triggers such a large difference in the 
estimate. 

We consider two sets of 10000 CMB sky simulations, 
each set generated using the same random seed, and 



f 


Tin 

7 


V+W 


V 


W 


Q 


324 


25 


64 (±71) 


77 (±73) 


63 (±75) 


26 (±74) 


390 


26 


44 (±61) 


81 (±64) 


35 (±66) 


25 (±66) 


471 


27 


44 (±55) 


71 (±60) 


40 (±62) 


31 (±62) 


567 


28 


42 (±52) 


56 (±56) 


55 (±58) 


43 (±56) 


683 


29 


73 (±49) 


71 (±52) 


72 (±54) 


23 (±50) 


823 


30 


81 (±43) 


80 (±46) 


72 (±49) 


24 (±46) 


1000 


31 


89 (±39) 


105 (±42) 


54 (±45) 


9 (±44) 



TABLE 2 

The estimated values for /jvl for different number rij of 

NEEDLET SCALES. SlNCE ERROR BARS INCREASE WHEN USING FEW 
NEEDLET SCALES, THE CORRESPONDING 1(T ERROR ESTIMATE IS 
GIVEN IN PARENTHESIS. 

therefore identical. One set is multiplied with the KQ75 
mask, while the other is multiplied with the KQ85 mask. 
Now we estimate /jvl for both sets, and find the dif- 
ference between each estimate, and the corresponding 
estimate from the identical map with the other mask, 

A/jvl = In? 75 ~ /jvj? 85 - Then the mean value and 
standard deviation of A/jvl is found. 

For 4 iax = 1000 we found Af NL = 28 for the WMAP 
data, whereas the standard deviation a&f NL = 21 for 
simulations. For £ max = 700, we found A/jvl = 68 and 



34. In the first case, the shift in /jvl when 



changing mask is as expected whereas in the latter case, 
the change A/jvl is slightly high, but only at the 2a 
level. 

As a further test of consistency, we also considered 
the difference in /jvl estimate between £ max = 700 



and £ ri 

A/iVL = / 



1000 when using the KQ85 galactic cut 



= 700 



JVL 



/ 



= 1000 



JVL 



= 39. However, a com- 
parison with simulations reveal that A/jvl have a stan- 
dard deviation of a — 30. In other words, the difference 
in the two estimates is well within 2a and is to be ex- 
pected. 

To test the variation of /jvl with increasing galactic 
cut, we estimated /jvl using the tiny Kpl2 mask as well 
as the extended KQ75± mask. We see that /jvl decrease 
when going from the smallest mask to KQ85 and KQ75, 
but increases slightly again to KQ75+. 

At this point a % 2 test was implemented. Three identi- 
cal sets of 10000 simulations were generated, and each set 
was multiplied with one of the KQ85, KQ75 and KQ75+ 
masks (we do not include the Kpl2 mask as foregrounds 
are likely to be important for such a small mask). For 
each simulation, a data vector, d, with two elements 
was formed from the difference in /jvl estimates between 
three different masks: 



f KQ75 _ rKQS5 
JNL J JVL 

pKQ75+ _ rKQ75- 
JNL JNL 



0) 



Of the 10000 simulations, 5000 were used to find a 
mean value and covariance matrix for d. Then a % 2 value 
was found for each of the remaining 5000 simulations as 
follows: 



X 2 = (d-(d)Y C- 1 (d-(d)) 



(10) 



A similar x 2 value was found from the /jvl values of the 
WMAP data maps. Then the x 2 values for the simula- 
tions were compared with that of the WMAP data. 

The result was that 37% of the simulations had a higher 
value of x 2 - We conclude that the variation of the /jvl es- 
timate for different masks (larger than Kpl2) are within 
expectations for a Gaussian map. 
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As a further check for possible foreground contamina- 
tion we will check the variation of Jnl with frequency 
channel. We investigated this by estimating /jvx using 
10000 simulated Gaussian CMB sky maps. For each sim- 
ulated sky, three identical maps were generated. These 
maps were then smoothed with the instrumental beam of 
the Q, V and W frequency channel respectively, and noise 
was added independently to each of the maps. For each 
of these sets of simulations a needlet transform was per- 
formed using 31 needlet scales in the range 2 < t < 1000. 
Then the bispectra were found and Jnl was estimated, 
using between 25 and 31 of the needlet scales. 

At this point we performed a % 2 test, similar to the one 
described above. First we tested the variation between 
the frequency channels when using all 31 of the needlet 
scales. For every simulation, a data vector, d, with two 
elements was formed from the difference in f^L estimates 
between the channels. 



L 



NL 



NL 



fV ' 
JNL 
fW 
JNL. 



(11) 



The results using the KQ75 mask showed that only 
4.5% of the simulations had a higher x 2 value, than the 
WMAP data. This corresponds to a « 2a deviation. 

To investigate whether this is consistent on several 
scales, we also performed the same test with some of the 
needlet scales removed. This was done using between 
25 and 31 needlet scales. Finally we combined the data 
vectors from all these tests: 



,Q31 
JNL 
fV31 
JNL ' 

f Q30 
JNL 
fV30 
JNL 



f Q25 
JNL 
fV25 
JNL ' 



fV31' 

' JNL 
fW31 
JNL 
fV3Q 

' JNL 
fW30 
JNL 



fV25 
' JNL 
fW25 
JNL ■ 



(12) 



and used this to make a combined test. From the KQ75 
results (shown in table [3]) it seems that only by using 
all available scales we find a small inconsistency of the 
Inl values between frequency channels. We repeated 
the latter test using the smaller KQ85 cut and found in 
this case that 14% of the simulations had a higher \ 2 
concluding that foreground residuals do not appear to 
be causing the difference in /jvl for different channels. 



To test the influence of foregrounds on the estimate of 
Jnl, we have estimated Jnl on the WMAP maps be- 
fore foreground subtraction (raw maps). The results are 
listed in table [TJ We see in particular for the Q band 
that the the value of /jvl is negatively biased by the 
pr esence of foregrounds. A similar result was also found 
in (jYadav & Wandelt ll2008HKomatsu et al.ll2008f) . Fore- 
ground residuals would thus be expected to give a too low 
value of Jnl- To check the power of our consistency test, 
we repeated the above x 2 test of the differences in esti- 
mated $nl between frequency channels using 31 scales. 
We find that only 0.7% of the simulations have a higher 
X 2 than for the WMAP data for the KQ75 cut, and none 
of the simulations have a similarly high x 2 for the KQ85 
cut. The test thus shows a clear detection of foreground 
residuals in this case. 

A similar x 2 test was now performed, but this time 
to study variation between different number of needlet 
scales (and thus also different i max ) used for the estima- 
tion: 



j-31 f30 ' 
JNL JNL 
f30 _ f29 
JNL JNL 



f26 

Jnl 



f25 
JNL. 



(13) 



where the superscript denotes number of needlet scales 
used to estimate the $nl value, x 2 was found using 
equation ((TOj) for the WMAP data as well as for the 5000 
simulations according to the same procedure as above. 
This was done for the individual Q, V and W frequency 
channels, and for the combined V+W map. The test 
was also performed using a combined data vector from 
all the three frequency channels. The results (table [4| 
show that the variation in the /n l estimate with respect 
to needlet scales is well within the expected bounds. 



Preq. channel 


% of sim. with higher \' 2 


Q 


91.9 


V 


44.3 


w 


72.2 


V + W 


76.6 


combined Q, V and W 


56.3 



TABLE 4 

Test of f NL variation with respect to scale. Fraction of 

SIMULATIONS WITH HIGHER \ 2 VALUE THAN WMAP DATA. THE 
RESULTS SHOW THAT THE WMAP DATA IS CONSISTENT WITH 
GAUSSIANITY IN THIS RESPECT. 





needlet scales 


% of sim. with higher \ 2 


324 


25 


20.6 


390 


26 


21.3 


471 


27 


54.6 


567 


28 


95.0 


683 


29 


43.2 


823 


30 


33.4 


992 


31 


4.5 


all of the above 


12.2 



TABLE 3 

Test of difference between frequency channels using 
different number of scales. the last row combines all 
the other variables in one test. the table shows 
percentage of simulations with higher \ 2 than the wmap 
data using the kq75 mask. 



5. CONCLUSIONS 

We have tested a n estimator for f^L bas ed on the 
needlet bispectrum (|Lan fc Marinuccil [20 08a). We used 
the estimator to obtain best fit values of Jnl from the 
WMAP 5 year data, using the combined V+W map as 
well as the independent frequency channels. The er- 
ror bars on Jnl obtained with the needlet bispectrum 
are significantly larger that those obtained by the opti- 
mal bispectrum estimator (|Smith et al. 1120081 ). but the 
needlet bispectrum still provides an important and in- 
dependent test of consistency. We have further intro- 
duced a set of consistency tests based on the difference 
A Jnl = /jvl — /jvx where 1 and 2 refer to different 
masks, different frequency channels or different number 
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of multipoles. We compare the differences A/atl for the 
different cases to the values obtained in simulations. 

We find our best estimate of l using the combined 
bispectrum from the V and W channels giving /jvl = 
76 ±38 using the KQ75 mask and £ max — 1000, consistent 
within la with the value of /jvl = 51 ± 32 obtained by 
the WMAP team as well as with the values ob tained by 
(lYadav k Wandelt I [20081 : iSmith et al. 1 f2008h all using 

4iax = 750. 

Using the combined V + W map and I max = 1000, we 
obtained f NL = 84 ± 40 for KQ75 and f NL = 110 ± 37 
using KQ85 (corrected for point source bias). This dif- 
ference in /nx using the two different masks was found 
to be within the 2a limit from simulations and thus con- 
sistent with expectations. In order to further limit the 
risk of foreground contamination, we estimated /nl on 
an extended KQ75 mask excluding 37% of the sky giving 
f NL = 97 ± 42. 

Using the independent frequency channels and the 
KQ75 cut, we obtained f^L = 9 ± 44 for the Q band, 
f NL = 105 ± 42 for the V band and f NL = 54 ± 45 for 
the W band. Such a large difference in f^L between the 
bands were found only in 4.5% of the simulations. This 
is 2(7 away from the expected value. This could be a sign 
of foreground residuals but could also well be a statistical 



fluke. We found the latter explanation to be more rea- 
sonable considering that for the smaller KQ85 mask 14% 
of the simulations had a larger difference. Similar tests 
were made with values of /nl obtained using different 
number of multipoles and channels, and no significant 
deviations from the expected differences were found. We 
therefore conclude that there are no convincing evidence 
of foreground residuals having influenced the estimated 
value of fNL, even using the KQ85 galactic cut. However 
repeating these tests on the next release of the WMAP 
data and on Planck data will be necessary in order to 
confirm this claim. 
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